3次元二値画像解析

このチュートリアルでは3次元の画像解析を行ないます。基本的な考え方は2Dの二値画像解析と同じなので、まずそちらのチュートリアルを終わらせてからこちらにトライしてください。

ここでは50x50x50のサイズの画像を解析します。3次元画像の表現方法は色々ありますが、ここでは各層ごとの断面画像を用います。 つまり50x50の画像を50個用意します。 data/ ディレクトリにファイルが置かれています。ファイル名は断面の順に付けておきます。

In [1]:
ls data/

0000.png  0007.png  0014.png  0021.png  0028.png  0035.png  0042.png  0049.png
0001.png  0008.png  0015.png  0022.png  0029.png  0036.png  0043.png
0002.png  0009.png  0016.png  0023.png  0030.png  0037.png  0044.png
0003.png  0010.png  0017.png  0024.png  0031.png  0038.png  0045.png
0004.png  0011.png  0018.png  0025.png  0032.png  0039.png  0046.png
0005.png  0012.png  0019.png  0026.png  0033.png  0040.png  0047.png
0006.png  0013.png  0020.png  0027.png  0034.png  0041.png  0048.png

データファイルの読み込み

まずはこの 50x50x50 の画像を読み込んで 50x50x50 の白黒配列を作ります。

In [2]:
# 必要なライブラリの読み込み
import imageio
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

3次元配列を作る前に 0000.png を読み込んでどんな画像か確認しましょう。

In [3]:
image_0000 = imageio.imread("data/0000.png")
In [4]:
image_0000.shape, image_0000.dtype, np.max(image_0000), np.min(image_0000)
Out[4]:
((50, 50), dtype('uint8'), 255, 0)

これは50x50の配列で、0 から255の値を持っています。つまりこの画像はグレイスケール画像のようです。中身を見てから実際に表示してみましょう。

In [5]:
image_0000
Out[5]:
Array([[255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       ...,
       [255, 255, 255, ...,   0,   0,   0],
       [255, 255, 255, ...,   0,   0,   0],
       [255, 255, 255, ...,   0,   0,   0]], dtype=uint8)
In [6]:
plt.imshow(image_0000, "gray")
Out[6]:
<matplotlib.image.AxesImage at 0x7f992e485130>

では50個の画像を読み込んで積み重ねます。

In [7]:
pict = np.stack([
    imageio.imread("data/{:04d}.png".format(n)) > 128
    for n in range(50) 
], axis=0)
In [8]:
pict.shape, pict.dtype
Out[8]:
((50, 50, 50), dtype('bool'))

50x50x50の画像が完成しました。ではこれを3次元可視化しましょう。plotly.graph_objectsと homcloud.plotly_3d モジュールを使います。 詳しくはAPIドキュメントを見てください。 Bitmap3d関数をここでは使います。コメントアウトの所を有効化すると不透明度を変えることができます．

In [9]:
# 必要なモジュールを読み込む
import homcloud.plotly_3d as p3d
import plotly.graph_objects as go
In [10]:
fig = go.Figure(data=[p3d.Bitmap3d(pict)], layout=dict(scene=p3d.SimpleScene()))
# fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.show()

パーシステント図の計算

パーシステント図を計算します。hc.PDList.from_bitmap_distance_functionを使います。 2Dの場合と同じです。

In [11]:
import homcloud.interface as hc
In [12]:
pdlist = hc.PDList.from_bitmap_levelset(hc.distance_transform(pict, signed=True), save_to="binary-3d.pdgm")

以下のようにしてファイルを読み込むこともできます。

In [13]:
pdlist = hc.PDList("binary-3d.pdgm")

パーシステント図の解析(1次)

1次のパーシステント図を調べていきましょう。

In [14]:
pd1 = pdlist.dth_diagram(1)
In [15]:
pd1.histogram((-15.5, 10.5), 26).plot(colorbar={"type": "log"})

(-4, 4) の所に2つbirth-death pairがあるようです。これを確認しましょう

In [16]:
pairs = pd1.pairs_in_rectangle(-4, -4, 4, 4)
pairs
Out[16]:
[Pair(-4.0, 4.0), Pair(-4.0, 4.0)]

確かに2つあります。

逆解析(1次)

この2つのペアがどのような穴を表現しているのか見てみましょう。1次のパーシステントホモロジーを見ているわけですから、 何らかのリング構造を見ているはずです。それを計算しましょう。

ここでは optimal cycle というものを使います。ポイントクラウドの解析ではoptimal volumeやvolume-optimal cycleというものを用いましたが、 それと似たものです(ちょっと違います)。Pair.optimal_1_cycleで計算します。

In [17]:
optimal_1_cycles = [pair.optimal_1_cycle() for pair in pairs]

まずはこれを可視化しましょう。元のボクセルデータを重ねあわせて表示します。

In [18]:
fig = go.Figure(data=[
    cycle.to_plotly3d() for cycle in optimal_1_cycles
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap"),
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.1, selector=dict(name="bitmap"))
fig.show()

回転させたりして観察してください．opacityを変更するなどしてみてみましょう． ボクセルデータの場合には1次のbirth-death pairに対応するのはトンネル構造である、というのが この可視化を見るとわかると思います。

次にこのリング構造がどういうパスを通っているのかを見ます。Optimal1CycleForBitmap.pathを使います。通過点のピクセルの座標が得られます。

In [19]:
optimal_1_cycles[0].path()
Out[19]:
[(49, 4, 44),
 (49, 4, 43),
 (49, 4, 42),
 (49, 4, 41),
 (49, 4, 40),
 (49, 4, 39),
 (49, 4, 38),
 (48, 4, 38),
 (48, 4, 37),
 (48, 4, 36),
 (47, 4, 36),
 (47, 4, 35),
 (47, 4, 34),
 (46, 4, 34),
 (45, 4, 34),
 (44, 4, 34),
 (43, 4, 34),
 (43, 4, 33),
 (42, 4, 33),
 (41, 4, 33),
 (40, 4, 33),
 (39, 4, 33),
 (38, 4, 33),
 (37, 4, 33),
 (36, 4, 33),
 (36, 5, 33),
 (35, 5, 33),
 (34, 5, 33),
 (34, 6, 33),
 (34, 7, 33),
 (34, 8, 33),
 (34, 9, 33),
 (34, 10, 33),
 (34, 11, 33),
 (34, 11, 34),
 (34, 11, 35),
 (34, 11, 36),
 (34, 11, 37),
 (34, 11, 38),
 (34, 11, 39),
 (34, 11, 40),
 (34, 11, 41),
 (34, 11, 42),
 (34, 11, 43),
 (34, 11, 44),
 (34, 11, 45),
 (35, 11, 45),
 (35, 11, 46),
 (36, 11, 46),
 (37, 11, 46),
 (37, 11, 47),
 (37, 11, 48),
 (38, 11, 48),
 (39, 11, 48),
 (40, 11, 48),
 (40, 10, 48),
 (41, 10, 48),
 (42, 10, 48),
 (43, 10, 48),
 (43, 9, 48),
 (44, 9, 48),
 (45, 9, 48),
 (45, 8, 48),
 (45, 8, 47),
 (46, 8, 47),
 (46, 7, 47),
 (46, 7, 46),
 (47, 7, 46),
 (47, 7, 45),
 (48, 7, 45),
 (48, 7, 44),
 (48, 6, 44),
 (48, 5, 44),
 (49, 5, 44),
 (49, 4, 44)]
In [20]:
optimal_1_cycles[1].path()
Out[20]:
[(16, 42, 45),
 (16, 41, 45),
 (16, 41, 46),
 (15, 41, 46),
 (15, 40, 46),
 (15, 40, 47),
 (14, 40, 47),
 (14, 39, 47),
 (14, 39, 48),
 (13, 39, 48),
 (13, 39, 49),
 (12, 39, 49),
 (11, 39, 49),
 (11, 38, 49),
 (11, 37, 49),
 (10, 37, 49),
 (10, 36, 49),
 (10, 35, 49),
 (10, 34, 49),
 (9, 34, 49),
 (9, 33, 49),
 (9, 32, 49),
 (8, 32, 49),
 (8, 31, 49),
 (8, 30, 49),
 (8, 29, 49),
 (8, 29, 48),
 (8, 29, 47),
 (8, 29, 46),
 (8, 29, 45),
 (8, 29, 44),
 (8, 29, 43),
 (8, 29, 42),
 (8, 29, 41),
 (9, 29, 41),
 (9, 29, 40),
 (9, 29, 39),
 (10, 29, 39),
 (10, 29, 38),
 (10, 29, 37),
 (10, 29, 36),
 (10, 29, 35),
 (10, 30, 35),
 (10, 31, 35),
 (10, 32, 35),
 (10, 33, 35),
 (11, 33, 35),
 (11, 34, 35),
 (12, 34, 35),
 (13, 34, 35),
 (13, 35, 35),
 (14, 35, 35),
 (14, 36, 35),
 (14, 37, 35),
 (14, 38, 35),
 (14, 39, 35),
 (14, 40, 35),
 (14, 41, 35),
 (14, 41, 36),
 (15, 41, 36),
 (15, 41, 37),
 (15, 42, 37),
 (15, 42, 38),
 (16, 42, 38),
 (16, 42, 39),
 (16, 42, 40),
 (17, 42, 40),
 (18, 42, 40),
 (18, 42, 41),
 (18, 42, 42),
 (18, 42, 43),
 (18, 42, 44),
 (18, 42, 45),
 (17, 42, 45),
 (16, 42, 45)]

パーシステント図の解析(0次)とその逆解析

次に0次のパーシステント図を解析していきましょう。

まず0次のパーシステント図をプロットします。

In [21]:
pd0 = pdlist.dth_diagram(0)
pd0.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})

(-16.0, -10.0), (-15.0, -10.0) という2のbirth-death pairを詳しく調べていきましょう。 ボクセルデータの0次と2次の逆解析にはBitmapPHTrees というクラスを使います(1次には使えないので、上で説明した optimal_1_cycle を代わりに使いました)。 BitmapPHTrees.for_bitmap_distance_functionはPDList.from_bitmap_distance_functionと似たインターフェースで 逆解析用の情報を計算します。

この PHTrees という名前は、ここで計算されるものが木構造であることを意味しています。 n次元のビットマップデータに対して0次とn-1次の2つの情報を同時に計算して、2つの木構造を持っています。 この木構造を取りだすためには PDList.bitmap_phtrees というメソッドを使います。

In [22]:
hc.BitmapPHTrees.for_bitmap_levelset(hc.distance_transform(pict, signed=True), 
                                     save_to="binary-3d-tree.pdgm")
Out[22]:
PDList(path=binary-3d-tree.pdgm)

以下のようにして計算したデータを読み込みます。

In [23]:
pdlist_with_tree = hc.PDList("binary-3d-tree.pdgm")

0次の情報を取り出すのは bitmap_phtrees(0) を使います。

In [24]:
phtree_0 = pdlist_with_tree.bitmap_phtrees(0)

(-16.0, -10.0), (-15.0, -10.0)という2つのbirth-death pairの情報を取り出すのには pair_nodes_in_rectangleを使いましょう。これは指定した範囲内にある birth-death pair (に対応する木構造のnode)をすべて取り出します。

In [25]:
nodes_0 = phtree_0.pair_nodes_in_rectangle(-16, -15, -10, -10)
nodes_0
Out[25]:
[BitmapPHTrees.Node(-15.0, -10.0), BitmapPHTrees.Node(-16.0, -10.0)]

取り出したデータは BitmapPHTrees.Node というクラスのインスタンスです。次にこれらの情報を可視化しましょう。

In [26]:
fig = go.Figure(data=[
    n.to_plotly3d() for n in nodes_0
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

0次のPHなので，ボクセルの内側の領域に 対応しています．実はこの2つのbirth-death pairは木構造上の親子関係にあり，一方が他方の部分集合になっています．よく見ると赤い領域の一部が青くなっているのが見えると思います．

パーシステント図の解析(2次)とその逆解析

最後に2次のパーシステント図を解析していきましょう。

2次のパーシステント図をプロットします。

In [27]:
pd2= pdlist.dth_diagram(2)
pd2.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})

(4, 6)というbirth-death pairを調べます。逆解析には bitmap_ph_trees(2)を使います。2が次数です．

In [28]:
phtree_2 = pdlist_with_tree.bitmap_phtrees(2)

(4, 6) という birth-death pairを調べましょう。まず pair_nodes_in_rectangle で (4, 6) という birth-death pair (に対応した tree の node) をすべて取り出します。

In [29]:
nodes_2 = phtree_2.pair_nodes_in_rectangle(4, 4, 6, 6)
nodes_2
Out[29]:
[BitmapPHTrees.Node(4.0, 6.0), BitmapPHTrees.Node(4.0, 6.0)]

2個の 以下が可視化した結果です。くりぬかれた領域の奥のほうがこのbirth-death pairに対応していることがわかると思います。

In [30]:
fig = go.Figure(data=[
    n.to_plotly3d() for n in nodes_2
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

volumeメソッドを使うことでこの領域の各座標を得ることもできます。

In [31]:
nodes_2[0].volume()
Out[31]:
[[27, 24, 25],
 [26, 24, 25],
 [27, 23, 23],
 [26, 24, 23],
 [27, 24, 23],
 [26, 23, 23],
 [26, 23, 24],
 [27, 24, 24],
 [26, 24, 24],
 [27, 23, 24],
 [27, 23, 25],
 [26, 24, 22],
 [27, 23, 26],
 [26, 24, 26],
 [27, 24, 26],
 [28, 25, 25],
 [28, 24, 25],
 [25, 25, 25],
 [26, 25, 25],
 [27, 25, 25],
 [28, 23, 25],
 [25, 23, 25],
 [26, 23, 25],
 [25, 24, 25],
 [27, 22, 25],
 [25, 22, 25],
 [26, 22, 25],
 [26, 25, 22],
 [27, 25, 22],
 [26, 22, 23],
 [25, 22, 23],
 [27, 22, 23],
 [25, 23, 23],
 [28, 24, 23],
 [25, 24, 23],
 [28, 23, 23],
 [28, 25, 23],
 [27, 25, 23],
 [26, 25, 23],
 [25, 25, 23],
 [27, 24, 21],
 [26, 24, 21],
 [27, 23, 21],
 [26, 23, 21],
 [27, 23, 22],
 [25, 24, 22],
 [26, 23, 22],
 [27, 24, 22],
 [28, 25, 24],
 [27, 25, 24],
 [26, 25, 24],
 [28, 23, 24],
 [25, 25, 24],
 [28, 24, 24],
 [25, 24, 24],
 [28, 22, 24],
 [25, 23, 24],
 [27, 22, 24],
 [26, 22, 24],
 [25, 22, 24],
 [29, 24, 25],
 [29, 25, 25],
 [24, 25, 25],
 [25, 26, 25],
 [26, 26, 25],
 [27, 26, 25],
 [28, 26, 25],
 [26, 24, 20],
 [26, 23, 20],
 [27, 23, 20],
 [27, 24, 20],
 [27, 22, 21],
 [28, 22, 21],
 [26, 22, 21],
 [25, 22, 21],
 [24, 25, 24],
 [29, 24, 24],
 [29, 25, 24],
 [25, 26, 24],
 [26, 26, 24],
 [27, 26, 24],
 [28, 26, 24],
 [28, 23, 21],
 [28, 25, 21],
 [25, 23, 21],
 [28, 24, 21],
 [25, 24, 21],
 [25, 25, 21],
 [26, 25, 21],
 [27, 25, 21],
 [29, 22, 24],
 [24, 23, 24],
 [27, 21, 24],
 [26, 21, 24],
 [25, 21, 24],
 [24, 22, 24],
 [28, 21, 24],
 [24, 24, 24],
 [26, 22, 19],
 [29, 23, 24],
 [26, 23, 19]]
In [ ]: